Freezing in random graph ferromagnets 
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Using T = Monte Carlo and simulated annealing simulation, we study the energy 
relaxation of ferromagnetic Islng and Potts models on random graphs. In addition to the 
expected exponential decay to a zero energy ground state, a range of connectivities for 
which there Is power law relaxation and freezing to a metastable state is found. For some 
connectivities this freezing persists even using simulated annealing to find the ground 
state. The freezing is caused by dynamic frustration in the graphs, and is a feature of the 
local search-nature of the Monte Carlo dynamics used. The Implications of the freezing on 
agent-based complex systems models are briefly considered. 

PACS numbers: 75.10.Hk [Classical spin models), 75.40.Mg (Numerical simulation 
studies), 75.10.Nr (Spin-glass and other random models) 
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The way a physical system approaches equilib- 
rium is a subject of interest to both physicists and 
mathematicians. In order to measure thermody- 
namlcal properties of systems it is important to 
be certain that the system really Is In equilibrium. 
To ensure this, in computer simulations using the 
Monte Carlo dynamics it is necessary to first run 
the simulation for a long time before measuring. 
The way that various properties (e.g., the energy) 
of the system change during equilibration is also 
interesting in itself, e.g., in studies of how an epi- 
demic disease or an opinion spreads in a model of 
social agents. 

Here we study the relaxation of the energy of 
ferromagnetic Ising and Potts models on random 
graphs using Monte Carlo simulations with the 
Metropolis dynamics. We find an interesting tran- 
sition as the connectivity of the graph is varied. 
For very small connectivities, the energy relaxes 
exponentially fast, for intermediate connectivities 
the system freezes in a local minimum (with power 
law relaxation to it), and for graphs with large con- 
nectivities there is again exponentially fast decay. 

The model under consideration is the stan- 
dard Ising model with ferromagnetic interactions 
but with spins placed on a random graph. In 
graph theory terminology |[l|] , the ensemble used is 
Q{N,M), which consists of all graphs with N ver- 
tices and M — ^jN randomly selected edges. On 
average, each node is connected to 7 others; 7 is 
the connectivity or average degree of the graph. 
Each edge in the graph is a ferromagnetic interac- 
tion between the two linked spins, and the energy 
of the model can be taken to be 
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where exactly M of the Jij 's are non-zero and equal 



to 1. Thus, e counts the number of edges linking 
spins with different values. Note that this differs 
from the standard ferromagnetic Hamiltonian by 
a 7-dependent term. Similar models on random 
graphs have been used to study many different 
systems in biology and social science as well as 
in physics (e.g., |g, |3[). 

The model can also be viewed as a constraint 
satisfaction problem. Each of the M edges in the 
graph corresponds to a constraint that the two 
linked spins should be equal. A natural interpre- 
tation of this is a model of a social system where 
there are N agents choosing from two different 
opinions or activities. A link between two agents 
would mean that the two prefer to agree. 

By relaxation of a model, we mean the behaviour 
of the energy after a quench from a high temper- 
ature disordered spin configuration. The Monte 
Carlo method tries to decrease the energy of 
the system by changing the configuration of spins 
locally. In the Glauber dynamics used in this 
paper, the change is accomplished by attempting 
to flip a randomly selected spin. If the new spin 
configuration has lower energy than the old, it is 
accepted. If the energy is raised A units by the 
change, the new configuration is accepted with 
probability cxp[— /3A], where /3 = 1/T is inverse tem- 
perature (this is the Metropolis Q algorithm). In 
temperature T = simulation, no changes that 
raise the energy are accepted. In most of the simu- 
lations reported here, the Mitchell-Moore additive 
generator (see, e.g., [^j was used to generate ran- 
dom numbers. Some runs were also performed us- 
ing the standard C library's drand4 8 ( ) generator; 
these gave the same results. 

For the standard 2D Ising model, with Jy = 1 if 
and only if spins i and j are nearest neighbours 
on the square lattice, two behaviours of the re- 




FIG. 1: The relaxation In a ferromagnetic Ising model 
on a random graph with 10"' vertices, averaged over 50 
graphs and 10 restarts per graph. For small (not shown) 
and large 7, there Is a fast exponential relaxation, whUe 
the behavior for 7 = 2 and 3 Is a power law e = eo +4^", 
with V PS 1.3. The lower arrow Indicates the curve for 
7 = 1.5 and the upper one the 7 = 8 data. In between 
them are data for 7 = 2, 3, 4, 5, 6, and 7. 



laxatlon are possible. If the order parameter is 
conserved by the dynamics, so that the magneti- 
sation of the system does not change, e t^^/^, 
while e ^ t^^/^ if single spin flip dynamics are used. 
These behaviours can be understood by consider- 
ing domains of up and down spins Q]. 

Since a random graph is locally tree-like, it is 
natural to approximate the behaviour of the ran- 
dom graph model with that of the same model on 
a tree. Johnston and Plechac |^] have shown that 
the thermodynamical behaviour of the ferromag- 
netic Ising model is independent of the presence 
of loops in a graph: it is the same on a random 
regular graph and on a Bethe tree with the same 
connectivity. Da Silva and Silva ||^] studied relax- 
ation in the Ising ferromagnet on Cayley trees, and 
showed that mean field theory predicts exponential 
relaxation. Exponential relaxation can also be ar- 
gued for easily by writing a mean-field equation for 
the time dependence of a spin in terms of its near- 
est neighbours. Glassiness in the Cayley tree fer- 
romagnet has been studied by Melin et al [ |To| ] , who 
find a crossover temperature that scales inversely 
with the logarithm of the number of surface sites. 
For the random graph model considered here, this 
is 0, since there are no surface sites. 

Figure H] shows the relaxation behaviour of e for 
7 = 1.5 to 8 and graphs of size 10"*. All data were av- 
eraged over 50 different graphs, and the MC algo- 
rithm was restarted in 10 different initial spin con- 
figurations for each graph. In order to check self- 
averaging, we also made runs with averages over 
5 graphs and 100 initial configurations, and 500 
graphs and 1 initial configuration, and found no 



differences. Error bars were determined to be on 
the order of 10^'' or smaller. The figure shows that 
large 7's cause very fast relaxation to the ground 
state, while the system freezes for intermediate 
values of 7 . For very small 7's, the relaxation is 
of course still fast (not shown in the figure). The 
behaviour for intermediate values of 7 is thus dif- 
ferent from the tree-like models. We have also ob- 
tained similar results using ferromagnetic k state 
Potts models on random graphs. 

This behavior can be explained qualitatively by 
noting that even though the ferromagnetic models 
always have a ground state with zero energy, it is 
possible for the T = Monte Carlo algorithm (and 
all other local search methods) to get stuck in a 
local minimum. The simplest case where this can 
happen is when there is a link between two nodes 
that have different values and each of the nodes 
have two neighbours with the same value, see fig- 
ure ^. Because there is only one path between the 
up and down domains in this figure, it is not pos- 
sible to lower the energy by fiipping a spin. Thus, 
even though the model itself is solvable and con- 
tains no frustration, the dynamics gives rise to dy- 
namical frustration for local search methods. (Very 
recently, Spirln et al |]_1] have found freezing to 
a blinker state in the 3-dimensional Ising model. 
Blinkers will appear in the random graphs studied 
here too, but because of the relatively small con- 
nectivities at which the freezing appears it is more 
likely that it is caused by subgraphs such as those 
shown in figure ^.) 

If there are sufficiently many edges between the 
up and down domains, the relaxation will be fast. 
No dynamical frustration will occur and one domi- 
nant value will spread quickly through the graph. 
If there are few edges, this will take longer, and 
different values will dominate different parts of the 
graph. This makes it plausible that introducing a 
metric and adding a restriction to the range of the 
edges could cause changes in the relaxation. This 
conjecture is confirmed by simulations of a model 
where the spins are arranged on a chain and edges 
only allowed between spins whose distance is less 
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FIG. 3: This figure shows e after 10'^ MC steps per spin 
as a function of 7 for system sizes ranging from 50 to 



FIG. 4: Some of the data from figure ^ rescaled by the 
maximum energy for each A'^ and plotted as a function 
of a rescaled parameter 7 described in the text. 



than aN, where a is Independent of N. The large 
7 behavior is now power law relaxation and freez- 
ing. This is similar to the behaviour of the anti- 
ferromagnetic Ising and Potts models on a random 
graph |12|. This is a model for the graph colour- 
ing problem, a combinatorial optimisation problem 
that is NP-complete [13], meaning that its worst 
case instances in all likelihood require exponential 
time to solve on a deterministic Turing machine. 

Returning to the model with no restrictions on 
the edges, figure ^ shows the value of e after lO'^ 
MC steps per spin as a function of 7 and for sys- 
tem sizes ranging from 50 to 10^, determined us- 
ing about 100 different graphs and about 50 runs 
for each graph. The freezing region can be seen 
clearly. The error bars of the results shown in this 
and the other figures were small, typically on the 
order of or smaller than the symbols used to plot 
the data. 

For large N, it is possible to fit all the data from 
figure ^ on a universal curve. Figure ^ shows the 
energy for large N, rescaled so that the maximum 
is 1 , as a function of a rescaled parameter 
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where 70 is the location of the maximum and 
was calculated so that 7 = ±1 marks the points 
where the energy attains half its maximum value. 
Note that the original data are non-symmetric 
around their maxima: when calculating 7 we di- 
vided by different factors right and left of the maxi- 
mum. To determine the locations of the maximum 
and i -maximum points, a cubic spline fit of the 
data was used; the plot shows the original data 
points. 

Since the energy barrier surrounding the local 
minimum shown in figure ^ is small, it is likely 
that finite temperature MC simulations would not 
show the same behaviour. To test this, we have 
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FIG. 5: This figure shows the energy after 10^ MC steps 
per spin for system sizes from 500 to 10'' using T = 
MC simulation (top left) and simulated annealing with 
start temperature Ti = 0.1, 0.2, and 0.4 (top right, bot- 
tom left, bottom right) and end temperature 0. Most of 
the freezing effect disappears using simulated anneal- 
ing, but even for the Ti — 0.4 runs the algorithm was 
unable to find the ground state for 7 2. 



also tried simulated annealing |[T^] on the problem. 
Simulated annealing starts at a high temperature 
and then gradually decreases it during the simu- 
lation. We used a linear decrease in temperature, 
T{t) ^Ti-kt, where k was chosen so that the sim- 
ulation ends at zero temperature. 

Figure |5| compares the values of e after 10'^ MC 
steps per spin for the T = MC algorithm with 
those obtained using simulated annealing with 
start temperature Ti — 0.1, 0.2, and 0.4 and the 
same averaging as in the previous runs. Most of 
the freezing disappears in these runs, but there 
is a region of 7's for which it remains. Note that 
the percolation threshold for random graphs is at 
7=1, well below the freezing. 

In conclusion, we presented results from Monte 
Carlo and simulated annealing studies of the fer- 
romagnetic Ising model on random graphs. We 
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find different regions of behaviour of e{t) — the 
expected exponential relaxation but also some re- 
gions where there is power law relaxation. More 
importantly, freezing was found in the model. The 
freezing persisted even for some simulated anneal- 
ing runs, but almost disappears for large start- 
temperatures. The freezing is a feature of the lo- 
cal search and hill-climbing characteristics of the 
MC method used. This has implications for the 
study of models of choice-making agents on ran- 
dom graphs: for some connectivities it is not pos- 
sible to reach a consensus or the globally most 



effective solution by using only local information. 
There are intriguing similarities and differences 
between this model and the corresponding antifer- 
romagnetic model studied elsewhere; these could 
be studied further by examining the model where 
there are mostly ferromagnetic bonds but with 
some probability p of instead having an antiferro- 
magnetic bond. 
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